!Calcul de \ln(N/N_now) \sim - \xi \ln(t)
!Calcul de \ln(v0/v) \sim - \gamma \ln(t)

module fonctions
implicit none
integer, parameter :: d=2 !dimension
integer, parameter :: N = 1000000 !nombre initial de particules
double precision, parameter :: p = 0.05d0 !probabilité d'annihilation. Attention: il faut p>0 (prévu: N->0). a_2(p=0)=0
double precision, parameter :: intervalle = 0.01d0 !intervalle log10(N/Nnow) de calcul des a_2
integer, parameter :: Nmin = 2 !nombre de particules minimales d'un run
integer, parameter :: Navg = 80000 !nombre de réalisations où toutes les particules sont annihilées pour moyenne
character*(*), parameter :: filename1 = 'a2ap-2xi.txt' !nom du fichier de sortie
character*(*), parameter :: filename2 = 'a2ap-2gamma.txt' !nom du fichier de sortie
double precision, parameter :: alpha = 1.0d0 !coefficient de restitution. Laisser \alpha=1.0
double precision, parameter :: pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164
contains

FUNCTION ran2(idum) !Initialize idum with negative integer. 1.5x +plus lent que ran1. ran2 si +de 100M nb aléatoires
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL :: ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,&
          &IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER :: idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if(idum.le.0)then
	idum=max(-idum,1)
	idum2=idum
	do j=NTAB+8,1,-1
		k=idum/IQ1
		idum=IA1*(idum-k*IQ1)-k*IR1
		if(idum.lt.0)idum=idum+IM1
		if(j.le.NTAB)iv(j)=idum
	enddo
	iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if(idum.lt.0)idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if(idum2.lt.0)idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
return
end function ran2

subroutine rescale(vitesses,Nnow) !enlève la vitesse du centre de masse et rescale la température
	integer, intent(in) :: Nnow
	double precision, dimension(1:Nnow,1:d), intent(inout) :: vitesses
	double precision, dimension(1:d) :: vcm !vitesse du centre de masse
	double precision :: dNnow,tmp
	integer :: l
	dNnow = dble(Nnow)
	vcm = 0.d0
	do l=1,Nnow
		vcm = vcm + vitesses(l,:)
	end do
	vcm = vcm/dNnow
	do l=1,Nnow
		vitesses(l,:) = vitesses(l,:) - vcm
	end do
!	tmp = 0.d0 !second moment = 1 (éviter que température -> 0)
!	do l=1,Nnow
!		tmp = tmp + sum(vitesses(l,:)**2)
!	end do
!	vitesses = vitesses*dsqrt(dNnow)/dsqrt(tmp)
return
end subroutine rescale

end module fonctions

!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------

program a2ap
use fonctions
implicit none
integer :: seed,i,j,k,l,p1,p2,Nnow,taille
double precision :: beta,norme,prodscal,omega,omegamax,xold,dN,dNnow,time,mu,v0
double precision, allocatable, dimension(:) :: x,t,g
integer, allocatable, dimension(:) :: nbmoyenne
double precision, dimension(1:N,1:d) :: v
double precision, dimension(1:d) :: sigma,vchange
logical :: collision
seed = -14731

dN = dble(N)
taille = int(dlog10(dN/dble(Nmin))/intervalle)+10 !première approximation
allocate(x(1:taille),g(1:taille),t(1:taille),nbmoyenne(1:taille))
x = 0.d0
g = 0.d0
t = 0.d0
nbmoyenne = 0 !tableau qui dit combien de valeurs j'ai pour moyenne du point i
beta = (1.d0 + alpha)/2.d0
avgmain: do l=1,Navg
	collision = .true.
	omegamax = 1.d0 !pas besoin de trouver un meilleur omegamax car les premières iterations ignorées l'améliorent
	k = 1
	Nnow = N
	dNnow = dble(Nnow)
	xold = 0.d0
	time = 0.d0
! ***************************************************************	
! *** Initialisation des vitesses entre 1 et -1 de façon isotrope
!	do i=1,N
!		distribisotrope: do
!			do j=1,d
!				v(i,j) = 2.d0*dble(ran2(seed))-1.d0
!			end do
!			if (sum(v(i,:)**2) < 1.d0) exit distribisotrope
!		end do distribisotrope
!	end do
! ***********************************************************
! *** Initialisation des vitesses en loi de Puissance en v^mu
! mu = -1.d0
!	do i=1,N
!		distribisotrope2: do
!			do j=1,d
!				v(i,j) = 2.d0*dble(ran2(seed))-1.d0
!			end do
!			if (sum(v(i,:)**2) < 1.d0) exit distribisotrope2
!		end do distribisotrope2
!	  do j=1,d
!	    v(i,j) = v(i,j)/(sum(v(i,:)**2)**(mu/(2.d0*(mu+dble(d)))))
!	  end do
!	end do
! ************************************************************************
! *** Initialisation des vitesses selon distrib. gaussienne ("Box Muller")
	do i=1,N
		distribgauss: do j=1,d
			v(i,j) = dsqrt(-2.d0*dlog(dble(ran2(seed))))*dcos(2.d0*pi*dble(ran2(seed)))
		end do distribgauss
	end do
!*************************************************************************
	call rescale(v(1:N,:),N)
	v0 = dsqrt(sum(v**2)/dN)

	boucle_principale: do
		!*** Créer au hasard un vecteur unitaire $\sigma$ de distribution isotrope
		normeiso: do !norme du vecteur doit être inférieure à 1 (rejette les vecteurs en dehors de l'hyper-cercle)
			do i=1,d
				sigma(i) = 2.d0*dble(ran2(seed))-1.d0
			end do
			norme = sum(sigma**2)
			if ((norme < 1.d0) .and. (norme > 0.d0)) exit normeiso
		end do normeiso
		sigma = sigma/dsqrt(norme)

		!*** Collisions
		differentparticles: do !choisit deux particules différentes au hasard
			p1 = int(ran2(seed)*dNnow + 1.)
			p2 = int(ran2(seed)*dNnow + 1.)
			if (p1 /= p2) exit differentparticles
		end do differentparticles
		prodscal = dot_product(v(p1,:)-v(p2,:),sigma)
		if (prodscal > 0.d0) then !collision seulement si les particules sont dans le même demi-espace
			omega = abs(prodscal) !\omega = ||v_{12} \cdot \sigma ||
			if (omega/omegamax > ran2(seed)) then !collision
				time = time + 1.d0/(dNnow**2.*omega) !TEMPS
				if (p > ran2(seed)) then !Annihilation: refait un tableau de vitesses v "sans trous" mais de dimension Nnow-2
					collision =.true.
					Nnow = Nnow - 2					
					dNnow = dble(Nnow)
					recopie: do !remet aux Nnow premières places les vitesses des particules: enlève les trous
						if ((p1 <= Nnow) .and. (p2 <= Nnow)) then
							v(p1,:) = v(Nnow+1,:)
							v(p2,:) = v(Nnow+2,:)
							exit recopie
						else
							if ((p1 > Nnow) .and. (p2 > Nnow)) then
								exit recopie
							else
								if ((p1 <= Nnow) .and. (p2 > Nnow)) then
									if (p2 == Nnow+1) then
										v(p1,:) = v(Nnow+2,:)
										exit recopie
									else
										v(p1,:) = v(Nnow+1,:)
										exit recopie
									end if
								else
									if (p1 == Nnow+1) then
										v(p2,:) = v(Nnow+2,:)
										exit recopie
									else
										v(p2,:) = v(Nnow+1,:)
										exit recopie
									end if
								end if
							end if
						end if
					end do recopie
					v(Nnow+1,:) = 0.d0
					v(Nnow+2,:) = 0.d0
				else !Collision (inélastique)
					vchange = beta*omega*sigma
					v(p1,:) = v(p1,:) - vchange
					v(p2,:) = v(p2,:) + vchange
				end if
			end if
			if (omega >= omegamax) omegamax = omega*1.1d0 !rescale le omegamax
		end if

		if ((dlog10(dN/dNnow)-xold >= intervalle) .and. (collision .eqv. .true.)) then !rép. log. valeurs
			collision = .false.
			call rescale(v(1:Nnow,:),Nnow) !important d'enlever VCM avant calcul de a2
			xold = dlog10(dN/dNnow)
			k = k+1
!*** Pour le calcul de \xi
			x(k) = x(k) + xold
!*** Pour le calcul de \gamma
			g(k) = g(k) + dlog10(v0/dsqrt(sum(v(1:Nnow,:)**2)/dNnow))
!*** Fin
			t(k) = t(k) + time
			nbmoyenne(k)=nbmoyenne(k)+1
		end if

		if (Nnow <= Nmin) then
			if ((k > taille) .or. (l==1)) taille = k	
			open(unit=1, file=filename1)
			open(unit=2, file=filename2)
			do i=2,taille !écriture de toutes les valeurs, premier terme t=0 à éliminer car on a \ln(t) = \ln(0)
				write(1,*) x(i)/dble(nbmoyenne(i)),dlog10(t(i)/(dble(nbmoyenne(i))))
				write(2,*) g(i)/dble(nbmoyenne(i)),dlog10(t(i)/(dble(nbmoyenne(i))))
			end do
			close(1)
			close(2)
			exit boucle_principale
		end if
	end do boucle_principale
	
	open(unit=1,file="state")
		write(1,'(" Run:",i5,"/",i5,", #p:",i4,", time:",f8.4,", x:",f8.4)') &
	& l,Navg,Nnow,time,x(k)/dble(nbmoyenne(k))
	close(1)
	v = 0.d0
end do avgmain
deallocate(x,t,g,nbmoyenne)
write(*,*)
write(*,*) "Programme termine."
read(*,*)
end program a2ap